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In this review paper we present a stable Lagrangian numerical method for com- 
^ puting plane curves evolution driven by the fourth order geometric equation. 

The numerical scheme and computational examples are presented. 

, Keywords: curve evolution, Willmore flow, surface diffusion, Lagrangian 

■ method 

in 

1. Introduction 

ON 

The main purpose of this contribution is to suggest a method for computing 
evolution of closed smooth plane curves driven by the normal velocity (3 
depending on the intrinsic Laplacian of the curvature k and curvature itself: 



(3 = -d 2 s k + b(k) (1) 

where b(k) is a C 2 function of the curvature k, b(0) = 0. Numerical aspects 
of evolution of plane curves satisfying (1) have been studied in 1 for the case 
of the surface diffusion flow with no lower order terms, i.e. f3 = —d^k, and 
in 2 for the case of the so-called Willmore flow for which the normal veloc- 
ity is given by (3 = — <9 2 fc — ifc 3 . Recall that the latter case corresponds 
to the motion of elastic curves, e.g. the model of Euler-Bernoulli elastic 
rod, which is an important problem in structural mechanics. The elastic 
curve evolution and surface diffusion can be found in many practical ap- 
plications as sintering (in brick production), formation of rock strata from 
sandy sediments, metal thin him growth etc. (see e.g. 3 ). 
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2. Governing equations 

We follow the so-called direct (or Lagrangian) approach. We represent a 
solution of (1) by the position vector x satisfying the geometric equation 
d t x = (3N + aT where N, T are the unit inward normal and tangent vec- 
tors. An immersed regular plane curve T can be parameterized by a smooth 
function x : S 1 — > M 2 , i.e. L = {x(u),u E for which g = \d u x\ > 0. 
Taking into account the periodic boundary conditions at u = 0, 1 we shall 
hereafter identify S 1 with the interval [0,1]. The unit arc-length parame- 
terization will be denoted by s, so ds = gdu. The tangent vector T and 
the signed curvature k of T satisfy T — d s x, k = det(d s x, d 2 x). Moreover, 
we choose the unit inward normal vector N such that det(T, N) = 1. Let a 
regular smooth initial curve To = Image(cco) be given. It turns out that a 
family of plane curves T t = Image(a;(. ,t)), t S [0, T), satisfying (1) can be 
represented by a solution to the following system of PDEs: 

d t k = -d^k + d 2 s b(k) + d s (ak) + k(k(3 - d„a), (2) 
d t r\ = -k/3 + d s a, g = cxp(^), (3) 

d t x = -dtx - k3 ~ k m d 2 s x + (a- '±d s (k 2 ))d s x , (4) 

subject to initial conditions fc(.,0) = k , g(.,0) = go, x(.,0) = x (.) and 
periodic boundary conditions at u = 0, 1 (cf. 4,5 ). 



3. Approximation scheme and numerical experiments 

The presence of a tangential velocity a in the position vector equation 
has no impact on the shape of evolving curves. As it was shown e.g. in 4-9 
for general curvature driven motions (nonlinear, anisotropic, with external 
forces) incorporation of a suitable tangential velocity into governing equa- 
tions stabilizes numerical computations significantly. It prevents the direct 
Lagrangian algorithm from its main drawbacks - the merging of numerical 
grid points and their order exchange. It also allows for larger time steps 
without loosing stability. In our numerical solution we consider tangential 
velocity given by a nonlocal tangential redistributions discussed in a detail 
m 4,5,8,9 j£ f u ows f rom 4 . 5 that the redistribution functional a satisfying 

d s a = kp - {kf3) T + u (L/g - 1) , (5) 

with a constant ui > 0, is capable of asymptotic uniform redistribution 
of grid points along the evolved curve. In our computational method a nu- 
merically evolved curve is represented by discrete plane points x\ where the 
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index i = 1, ...,n, denotes space discretization and the index j = 0, ...,m, 
denotes a discrete time stepping. The linear approximation of an evolving 
curve in the j-th discrete time step is thus given by a polygon with vertices 
x J iy i — 1, ...,n. Due to periodicity conditions we shall also use additional 
values tLi — a^_i, x$ — x J n , x 1 n+1 = x\, x J n+2 = x J 2 . If we take a uniform 
division of the time interval [0, T] with a time step r = ^ and a uniform di- 
vision of the fixed parameterization interval [0, 1] with a step h = ^, a point 
a:^ corresponds to x(ih,jr). The systems of difference equations correspond- 
ing to (2)-(4) and (5) will be given for discrete quantities a-, rjf, rj, kj, 
xj, i = 1, ...,n, j = 1, ...,m, representing approximations of the unknowns 
a, n, gh, k, and x, respectively. Here aj represents tangential velocity of 
a flowing node x\, and rjf, rj w — kf, v\ represent piecewise 

constant approximations of the corresponding quantities in the so-called 
flowing finite volume x\_ x ,x\ . We shall use the corresponding flowing 

dual volumes x{_ ll SP i , where x\ = x *- 1 ^~ x * ; w ith approximate lengths 

qf w |^ — At the j-th discrete time step, we first find discrete values 

of the tangential velocity a? by discretizing equation (5). Then the values of 
redistribution parameter ijf are computed and utilized for updating discrete 
local lengths r\ by discretizing equations (3) . Using already computed local 
lengths, the intrinsic derivatives are approximated in (2), and (4), and pen- 
tadiagonal systems with periodic boundary conditions are constructed and 
solved for discrete curvatures kj and position vectors xj. In the sequel, we 
present in a more detail our discretization. Using r\~ x as an approximation 
of the length of the flowing finite volume x{Zitx{~ 1 at the previous jth 
time step we construct difference approximation of the intrinsic derivative 
d s a « — j^T 1 and by taking all further quantities in (5) from the previous 
time step. We obtain the following expression for discrete values of the tan- 
gential velocity: a{ = a{_i + r\~ k\ 1 0\~ — r\ 1 B^ 1 + ^(M- 5-1 — r\ ) 

where ^ = --^ {^S^ - !^^-Cj+ b ^),^ = 

n n 

Mi- 1 = ±Li-\ V- 1 = <£tj- \ B?- 1 = j+rr E^" fc r A J " and 

1=1 1=1 
oPq = 0, i.e. the point x 3 is moved in the normal direction only. In- 
serting (5) in (3) and using a similar strategy give us: r\ 1 — — ^ — = 
-r?~ .B- 7-1 + w(M J_1 - rj^ 1 ), for i = l,...,n,rf = rf n , rf n+1 = n{. 
Next we update local lengths by the rule: rj — exp(?7^),rl 1 = r J n _ 1 , Tq = 
r in r n+i = r i> r n+2 = r 2 • Subsequently, new local lengths are used 
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for approximation of intrinsic derivatives in (2) and (4). First, we derive 
a discrete analogy of the curvature equation (2). We have to approxi- 
mate the 4-th order derivative of curvature inside the flowing finite volume 
[xi-\,Xi], i = 1, ...,n. For that goal we take the following finite difference 
approximation: dtk(xi) rts k i+2 — ( 1 r h -tt + 

-5-^ )k i+ l + ( it r-^ + T-^ h-s4 1 2^ )fci+(-5-^ h 

tt 1 1 )ki—t H ki-2- Approximat- 

rfqf_ 1 r i qf_ 1 r i - 1 r i q i - 1 r i - 1 q i -2 > 1 1 r i g i _ir i _ig i _ 2 1 z 11 

ing first and second order terms in (2) by central differences and taking 
semi-implicit time stepping we obtain following pentadiagonal system with 
periodic boundary conditions for new discrete values of curvature: 



4H-2 + HH-i + 4H + •':><. , + e{k{ +2 = a 

subject to periodic b.c. k 3 _ x = k :l n _ 1 ,k 3 — k 3 n , k° n+1 = k{ 7 k J n+2 — k 2 , where 



9i-i r i-i9i-2 9i r i+i1i+i 

fj = + HH+l) - KM' 1 ) KM' 1 ) KCi ) 



T <A fli-i 

\P. = - I —. V — ; 1 ; : 1- — : — ) + -5— - 



,112 1 1 

C ■ = : : 1 : : 1 : : : 1 : : 1 : : h 

(4) 2r i+i r i(4) 2 r iQi1i-i rf(9i_i) 2 (Qf-i)M-i 
r r< *< ^ + 2 2 ' 



In order to construct discretization of equation (4) we approximate the 
intrinsic derivatives in a dual volume For approximation of the 

fourth order intrinsic derivative of the position vector we take similar ap- 
proach as above for curvature, but in the middle point x; t of the dual volume. 
In such a way and using the semi-implicit approach, we end up with two 
tridiagonal systems for updating the discrete position vector: 



A3 rJ I nj ,,,3 i no „j I ivj I cj j — -cj 
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for i = 1, ...,n subject to periodic b.c. x J _ 1 = x 1 n _ 1 , x 3 — x J n , x 3 n+1 = 

x l> X n+2 = x 2 where 

1 1 

A 3 — f 3 - 



7 7 1 1 1 7 7 7 ' 

' i-1 ' i+W+V i+2 

B 3 = -(^^ + -^ + ' + '\ + 



r 



i+1 

m)+^H-i) 4 3 (kj +1 ) 2 -(M) 2 

2r\ 2 4 

^7/1 1 1 1 \ 

V\ = - — 1 : H : , h — , : : + 

\ r kM+i ( r l+i) 2 4 (^+i) 2 9i+i r l+i4+i r l+2 ) 



2r? +1 2 4 ^ 

= ^ - (Ai+Bi+vi+si), n = q M~\ 

T T 



where </>(fc) = k 2 — -^-. The initial quantities for the algorithm are computed 
from a discrete representation of the initial curve xq, for details see. 5 Every 
pentadiagonal system is solved by mean of Gauss-Seidel iterations. Next 
we present results of numerical simulations for the curve evolution driven 
by (1). In our experiments evolving curves are represented by n = 100 
grid points and we use discrete time step r = 0.001. First wc numerically 
compute time evolution of the initial ellipse with the halfaxes ratio 2:1 for 
the case b(k) = (surface diffusion flow). We consider the time interval 
[0, 2]. The evolution of curves without considering tangential redistribution 
indicates accumulation of some curve representing grid points and a poor 
resolution in other parts of the asymptotic shape, see also Figure 1, a). 
In the case of asymptotically uniform tangential redistribution, we can see 
a uniform discrete resolution of the asymptotic shape, see Figure 1, b). 
Next we present evolution of an nontrivial initial curve driven by (1) with 
b(k) = 0. We show evolution of a highly nonconvex initial curve (see Figure 
1, c) with asymptotically uniform redistribution (oj = 1). Since elastic curve 
dynamics is very fast in case of highly varying curvature along the curve we 
have chosen smaller time step r = 10~ 6 . In Figure 2 we present evolution 
of an initial asteroid driven by (1) with b(k) = —^k 3 (the Willmore flow). 
A solution computed by the direct Lagrangian method is depicted by cross 
marks whereas solid curves correspond to the solution computed by the 
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level set method approach. For details we refer the reader to. 2 




a) b) c) 

Fig. 1. Surface diffusion of an initial ellipse and its circular asymptotic shapes without 
a) and with b) tangential redistribution. Evolution of the highly nonconvex initial curve 
and its asymptotic circular shape at time t = 0.170 c). 




Fig. 2. An asteroid as an initial condition and its evolution by the Willmorc flow at 
times t = 0, 0.0005, 0.005. 
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